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Abstract 

We present the study of ten random realizations of a density field character- 
ized by a cosmological power spectrum P{k) at redshift z = 50. The reliability 
of such initial conditions for iV-body simulations are tested with respect to their 
correlation properties. The power spectrum P(k), and the mass variance ^(r) 
do not show detectable deviations from the desired behavior in the intermedi- 
ate range of scales between the mean interparticle distance and the simulation 
volume. The estimator for £(r) is too noisy to detect any reliable signal at the 
initial redshift z = 50. The particle distributions are then evolved forward until 
z = 0. This allows us to explore the cosmic variance stemming from the random 
nature of the initial conditions. With cosmic variance we mean the fact that 
a simulation represents a single realization of the stochastic initial conditions 
whereas the real Universe contains many realizations of regions of the size of the 
box; this problem affects most importantly the scales at about the fundamental 
mode. We study morphological descriptors of the matter distribution such as 
the genus, as well as the internal properties of the largest object(s) forming in 
the box. We find that the scatter is at least comparable to the scatter in the 
fundamental mode. 
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1 Introduction 

Our present understanding of the formation and properties of the cosmological large- 
scale structure relies to a large extent on A-body simulations: given the difficulty 
in addressing theoretically the highly nonlinear regime of the growth of density inho- 
mogeneities by the gravitational instability, simulations have proven a valuable tool 
to get insight into the (non-linear) structure formation scenarios. Therefore, it is of 
considerable importance to confirm the reliability of such simulations. 
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It has been claimed recently (Baertschiger & Sylos Labini 2002) that there are 
major problems with generating initial conditions (ICs) for iV-body simulations. We 
can identify several reasons why the ICs may introduce uncertainties in the subsequent 
evolution. First, there is the problem of finite-mass resolution or discreteness: the 
initial continuum density field is modeled by the distribution of a finite number of 
point particles N, therefore only a finite number of Fourier modes of the density field 
can be reproduced reliably. The maximum wavenumber (Nyquist wavenumber) is given 
by &max = vr/Ax, where Ax is the mean interparticle separation. The modes k ^ k max 
have spurious values related to the point-particle distribution and may lead to artificial 
effects in the posterior dynamical evolution. The finite-mass resolution is expected to 
be irrelevant if the nonlinear mode-mode coupling to the modes k ^ k max has only a 
small influence on the dynamics. 

The second problem with the ICs is the finite size of the simulation box with side 
length B, which implies that the values of the Fourier modes of the density field with 
wavenumber smaller than the fundamental wavenumber, k < k m[n = 2n/B, are artifi- 
cially set to zero. This leads to two possible difficulties: first, the absence of mode-mode 
coupling to those large-scale modes, and second the so-called cosmic variance, meaning 
that the simulation box represents only one (finite-sized) realization of the stochastic 
initial density field, whereas the true Universe contains many realizations of regions of 
the size of the box. Therefore, the morphological properties of the matter distribution 
in a certain volume, as measured by e.g. the genus statistics, will presumably show 
some intrinsic scatter when placing the volume at different locations in the real Uni- 
verse. And this will also happen with the (internal) properties of any given class of 
objects. This is one of the main aspects of the current study and what we refer to as 
cosmic variance (in iV-body simulations) throughout the paper even though one might 
argue that this is not the "real" cosmic variance but rather an artificially introduced 
sampling variance. However, we are actually interested in exactly that (sampling) ef- 
fect which can easily be tested by simply running the same cosmological simulation 
but using different random realizations of the initial density field. 

In this work we study systematically the reliability of the initial density field used 
as an input to the iV-body simulations as well as the effect of their random nature onto 
the internal properties of clusters. In Section 2 we briefly review the most commonly 
way used to generate ICs and the code used to evolve the particles into the non- 
linear regime. Section 3 focuses on some of the statistical characteristics of the Dark 
Matter field: we consider the 2-point correlation function, the power spectrum, the 
mass variance in spheres, and the Minkowski functionals, the latter being sensitive to 
correlations of higher order. Finally, in Section 4 we investigate dark matter clusters 
identified within the simulations and quantify the effect of the cosmic variance on their 
internal properties. 
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2 TV-body Simulations 



2.1 Generating Initial Conditions 

The commonly used way for setting up initial conditions for a cosmological simulation 
is to make use of the Zeldovich approximation to move particles from a Lagrangian 
point q to an Eulerian point x (e.g. Efstathiou, Frenk & White 1985): 

x = q-D(t)S(q) , (1) 

where D(t) describes the growing mode of linear fluctuations and S(q) is the 'dis- 
placement field'. This method is not restricted to a cosmological scenario nor to the 
Zeldovich approximation: it is very general, relying only on the continuity equation for 
the transport of particles in the limit D(t) — > 0. The initial Lagrangian coordinates q 
are usually chosen to form a regular, three-dimensional lattice although there are other 
possible point-particle realizations yielding a homogeneous and isotropic density field 
on large scales (i.e. glass-like initial conditions, White 1996). 

For the runs presented in this study we used the code described in Klypin & Holtz- 
mann (1997) to set up the initial conditions 

S(q)=V q $(q), $(<f) = $> jg cos(fc-$ + & jg sin(fc-$ , (2) 

k 

where the Fourier coefficients and are related to a pre-calculated input power 
spectrum of density fluctuations, P(k), as follows: 

a S = Ri^y/m, b^^y/m- (3) 

Ri, i?2 are (Gaussian) random numbers with mean zero and dispersion unity. The 
factor 1/k 2 is (the Fourier transform of) the Green's function of Poisson's equation 1 and 
$(<f) can therefore be understood as the gravitational potential created by a Gaussian 
stochastic density field whose power spectrum agrees with the input P(k); the power 
spectrum P(k) measures the strength of each individual /c-mode contributing to the 
density field. However, to fully preserve the random nature both amplitudes (sine- and 
cosine-wave) are to be picked from a Gaussian distribution. 
Eq. (2) can be rewritten introducing complex numbers: 

= Y, A k e M*[k-q + dk\), ^exp(^) := ±[a s + a_ n - i{b u - b_ % )\. (4) 

k 

Both Aj: and % need to be drawn from appropriate random distributions. However, 
the ICs of cosmological relevance are ergodic for A^, with k ^> k min , making their 
random nature irrelevant: spatial regions of size much smaller than the simulation 

1 Actually, —1/k 2 is the correct Green's function, but the factor —1 can be dropped as Ri and i?2 
scatter around zero. 
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box already work as many different realizations inside the box for those amplitudes. 
Thus, cosmic variance enters through the random nature of the phases and of the 
amplitudes for k ~ k min . 

The idea of this paper is to create a certain number of random realizations of the 
same power spectrum P(k) by using different random seeds when drawing Ri, R 2 in 
Eq. (3). The input power spectrum P(k) was calculated using the CMBFAST code 
(Seljak & Zaldarriaga 1996), and all parameters (e.g. box size, number of particles, 
force resolution, integration steps, etc.) were fixed except for the seed for generating 
the random sequence providing the R values 2 . 

2.2 Simulation Details 

We created a data base of ten simulations that all were started at a redshift z% = 50 and 
evolved until z = in a ACDM (fi = 0.3, fi A = 0.7, Q 6 /i 2 = 0.04, h = 0.7, a 8 = 0.9) 
cosmological model using 128 3 particles within a box of side length B = 64/i -1 Mpc, 
giving a mass resolution of m p = 1.04 • lO 10 /*" 1 M . They were carried out using 
the publicly available adaptive mesh refinement code MLAPM (Knebe, Green & Binney 
2001). MLAPM reaches high force resolution by refining all high-density regions with 
an automated refinement algorithm. The refinements are recursive: the refined regions 
can also be refined, each subsequent refinement having cells that are half the size of the 
cells in the previous level. This creates an hierarchy of refinement meshes of different 
resolutions covering regions of interest. The refinement is done cell- by-cell (individual 
cells can be refined or de-refined) and meshes are not constrained to have a rectangular 
(or any other) shape. The criterion for (de-)refining a cell is simply the number of 
particles within that cell and a detailed study of the appropriate choice for this number 
can be found elsewhere (Knebe, Green & Binney, 2001). The code also uses multiple 
time steps on different refinement levels where the time step for each refinement level 
is two times smaller than the step on the previous level. A regular 256 3 domain grid 
was used to cover the whole computational volume in all runs and cells were refined as 
soon as the number of particles per cell exceeded the preselected value of 8. We stored 
snapshots of the particle distribution at redshifts z — 5, z — 1, z — 0.5, z = 0.25, 
and z = 0. At the end of the runs the force resolution is determined by the highest 
refinement level reached: for the runs at hand the finest grid at z = consisted of 8192 
cells per side and was called into existence at redshift z ~ 0.88. This grid corresponds 
to a force resolution of about 23h~ 1 kpc which is simply three times the grid spacing 
and gives the scale where the forces are purely Newtonian. This is sufficient enough for 
the presented study as we are mainly interested in the overall (large-scale) clustering 
properties. But as we will see later on, we are resolving approximately 2% of the virial 
radius of the most massive halo formed in the runs, which is sufficient for investigations 
of the internal properties such as velocity dispersion, spin parameter and triaxiality. 

2 An appropriate routine might be GASDEV from Numerical Recipes (Press et al. 1992). 
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Figure 1: Power spectrum evolution for all ten runs as compared to the prediction by 
Peacock & Dodds (1996) (solid line) and the linear P(k) (dotted line), respectively. 
The inset panel for z = 50 focuses on the fundamental mode k min = 2ir/B which shows 
a la variance of approximately 20%. 

3 Analysis I: Dark Matter Field 

We first focus on the properties of the dark matter particle distributions. Our main 
aim is to assess the recent claims by Baertschiger & Sylos Labini (2002) that there 
are major problems with generating initial conditions for iV-body simulations in the 
way as outlined in Section 2.1. Because it is common to either use a regular grid or a 
glass-like distribution as Lagrangian starting points q for the Zeldovich approximation 
(cf. Eq. (1)), their arguments try to prove that this leads to spurious artifacts related 
to e.g. the regular structure of such a grid, and that the initial conditions are not able 
to reflect the superposed CDM-like fluctuations at all. 

3.1 Power Spectra 

When creating a fluctuating density field in a certain volume by using a fixed number 
of particles, one is limited in the range of fc's by the size of that volume on the one 
hand, and the number of particles used to sample the waves on the other hand. The 
wave number of the lowest frequency wave (fundamental mode) to fit into the box is 
given by k min = 2ir/B where B is the side length of the box. The maximum wave 
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number is determined by the Nyquist frequency, k max = tt/Ax, where Ax = B/N 1 ^ 
is the mean particle separation (not to be confused with the grid spacing used in the 
iV-body code and for extracting the power spectrum from such a particle distribution, 
respectively). A recent investigation showed that high-resolution iV-body simulations 
where even smaller scales than k max are resolved are justified for power spectra with an 
effective spectral index n eff = d log P(k)/d log k much less than -1 (Hamana, Yoshida & 
Suto 2001). And this is the case for (nearly) all CDM type spectra as P(k) oc k~ 3 for 
large k. The evolution of power on small scales is driven by the transfer of power from 
large scales and hence it is important to follow that evolution with an appropriate force 
resolution even though that small-scale power was not present in the initial conditions 
(see Introduction). 

Using the particle data at redshifts z = 50 (initial conditions), z — 5, z — 1, and 
z = 0, we derived P(k) by Fourier transforming the density field on a regular 512 3 
grid 3 , which effectively introduces k max ~ 25h Mpc -1 as the maximum wavenumber 
to be recovered from the data. We adopted the method for extracting even higher k 
waves from the particle distributions as outlined in Jenkins et al. (1998). The power 
spectra were then compared to the non-linear prediction given by Peacock & Dodds 
(1996, PD96). However, their fitting parameters depend on the spectral slope n = 
d\nP(k)/dlnk and hence some recipe needs to be adopted when applying it to a 
cosmological P(k) where n is a function of k. We used n e s defined via 



d\nP(k) 

n e s{ki) = 



dink 



(5) 

k=k l /2 



for the estimate of the spectral index n cS at wave number k t (cf. Jenkins et al. 1998 
and Peacock & Dodds 1996). 

In Fig. 1 the results are shown along with the linearly extrapolated P{k). There 
are a couple of things to note besides the overall good agreement of the estimated 
P(k) with the PD96 prediction: first, the power spectrum derived from the particle 
distribution agrees extremely well with the input P(/c) 4 and the fundamental mode has 
not turned non-linear at z = 0; second, we can clearly see how the scatter prominent in 
the large waves k min ^ 2n/B at high redshifts moves towards higher k- values at later 
times. The scatter in the initial conditions is of the order of 20% and it arises because 
only a small number of such harmonics do fit into the finite box of side length B. We 
ascribe the migration of the scatter downwards to smaller scales (higher fc-values) to 
the transfer of power from large to small scales: the higher the amplitude at k min 
(cf. Eq. (4)) the more power can be transfered to smaller fc's. And hence we are facing 
a faster evolution of small-scale structures leading to the observed dispersion amongst 
the individual runs. 

The discrepancy with the PD96 prediction for z = 5 for k > 10 is not physical. 
Even though the MLAPM code invoked already three levels of refinements at z = 5 they 
are still very small in size, i.e. there are only ca. 40,000 refinement cells in total with 

3 The density was assigned to the grid cells using the Triangular-Shape-Cloud method. 

4 At redshift z — 50 the non- linear P(k) as given by PD96 is indistinguishable from the linear P{k). 
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about 12,000 particles (~ 0.6% of all particles) being moved on those levels. However, a 
visual inspection of the refined regions shows that the grids are covering all prospective 
halo formations sites and hence we are following the built-up of structures correctly. 
But when trying to recover those high-/c modes from the simulations we are left with 
observed mismatch due to the majority of the particles still being moved on the 256 3 
grid. 



3.2 Mass Variance a(r) 

The variance a M of the mass in spheres with radius r is given by 

a 2 M (r) = — / P(k)W 2 (kr)k 2 dk, W(x) = — (sin x - x cos x) . (6) 

2,71 Jo X 6 

The function cr| f (r) is readily calculated and can be compared to an adequate estimator 
<T Mest( r ) when being applied to the actual particle data. Our estimator distributes a 
certain number of spheres with radius r at random in the simulation volume and 
compares the number of particles inside those spheres to the expected mean value: 

2 ( \ 1 f (N t (r) - (N r )) 2 

^M,estW - ^ — - • [7) 

N s is the total number of spheres with radius r and (N r ) = (p)47rr 3 /3m p is the mean 
number of particles in such a sphere. 



3.2.1 Reliability of Estimator 

In order to make sure our estimator works as expected we started by applying it to 
particle distributions where simple scaling laws for (J 2 M {r) can be calculated analytically. 
For a purely Poissonian particle distribution one derives easily 

^M.PoissonM OC r" 3 , (8) 

and for a "shuffled" lattice (e.g. Gabrielli, Joyce & Sylos Labini 2002): 

^Lattice ( r ) oc r" 4 , (r > lattice spacing) . (9) 

From Fig. 2 we deduce that our estimator does indeed work correctly: we created 
ten Poisson distributions of 128 3 particles in a (128/?," 1 Mpc) 3 volume and for each 
distribution we calculated cr|^ cst (r) using 10000 spheres (for each r value). The curve 
shown is the average taken over the ten Poisson distributions. The error bars are 
too small to be presented. The shuffled lattice distribution was created as follows: we 
placed 128 3 particles on the nodes of a 128 3 grid with spacing l/i _1 Mpc, and then each 
particle was shifted in x, y, and z directions by a random amount uniformly distributed 
in the range [—0.05, 0.05] h~ x Mpc; ten such realizations were created. The curve shown 
in Fig. 2 is again the mean estimate when averaging over the ten realizations. In both 
tests we recover the expected scaling relation. 
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Figure 2: Reliability check for our o 2 Mcst (r) estimator Eq. (7). The solid lines have 
the slopes of the analytical expectations (refer to the text for further details). All 
amplitudes are arbitrary. 

3.2.2 Application to A-body Data 

We now apply the estimator Eq. (7) to our ICs as well as the final outputs at redshift 
z — 0. Fig. 3 shows j^^r) compared to the analytical cr^(r) as given by Eq. (6). 
For every scale r we used again N s = 10000 randomly placed spheres. The mean mass 
variance (o"M iCSt ( r ))set (averaged over the ten realizations) is plotted and the error bars 
are 1 times the variance of cr 2 MfiSt {r-) around the mean value (o"M ie st( r ))set- 

Contrary to the findings of Baertschiger & Sylos Labini (2002), we observe that 
the initial conditions agree from approximately the scale of the particle Nyquist fre- 
quency out to nearly half the box size with the analytical predictions. The faster drop 
°f (°M,est( r ))set for scales approaching the box size is simply the effect of the finite 
(periodical) box. As soon as the volume of the sphere comes close to the actual box 
size (which happens for r ps B/2) one finds nearly all particles in the sphere due to the 
periodic boundary conditions. Hence the variance a 2 M est (r) drops faster than predicted 
by Eq. (6). And the larger amplitude of o 2 M ^{f) for small scales is indeed a reflection 
of the discreteness of the initial conditions. But in any case Fig. 3 is a rather convincing 
argument that the mass variance in the initial conditions does agree with the CDM 
type fluctuations as described by the input power spectrum P(k). 

When comparing o"^ est (r) for the final output at redshift z = to the analytical 
ovf/-(r) in Fig. 4 (both using the linearly extrapolated P(k) as well as the non-linear 
P(k) given by PD96) we notice again a couple of things: firstly, the large scales (r ^ 
B/2) are still below the expectation, and secondly, there is more pronounced scatter 
on scales r < O.Qh^ 1 Mpc than found for the ICs. For r ^ B/2, the explanation 
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Figure 3: Mean value of cr£f(r) as given by Eq. (7) when averaged over the ten initial 
conditions at redshift z = 50. The error bars are la errors. The solid line is the ana- 
lytical expectation given by Eq. (6). The vertical line indicates the scale corresponding 
to the (particle) Nyquist frequency, k max . 




Figure 4: Same as Fig. 3 but this time for the final output at z — 0. The two analytical 
curves are the linear extrapolation of Eq. (6) to z = (dotted line) and the prediction 
for a 2 (r) when using the PD96 power spectrum with Eq. (6) (solid line). 
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is again the finite periodic box. The increased value for the variance cr| f est (r) for 
small scales r < 1/tT 1 Mpc (and its large scatter) is readily explained by the fact that 
gravitationally bound objects (and voids) are forming which introduces some sort of 
"semi-discreteness": this gives rise to a higher variance (as well as larger scatter) on 
scales related to the average size of such objects, i.e. Mpc and below. 



3.3 Two-Point Correlation Function £(r) 

The two-point correlation function is the Fourier transform of the power spectrum: 

« r '=ir pw ^ 2 *- <io) 

The basic interpretation of £(r) is that it is the average number of neighbors to a given 
object with distance r in excess to a Poisson distribution. And this is how we realize 
an estimator for £(r). We start again putting down a certain number of spheres in the 
simulation box, but this time centered at particles. We then create a shell of thickness 
dr extending from r to r + dr. The correlation function can now be estimated by 

U (r) = ^ - 1 , (ID 

where (p) is the mean number density of the simulation and T(r, dr) the mean number 
density of particles found in the shell [r, r + dr]: 

I n 3 

T(r,dr) = —J2Ur,dr) . (12) 

For each value of r, we use again N s = 10000 shells [r, r + dr] centered at a randomly 
chosen particle for calculating the average number density of particles T(r, dr) . 



3.3.1 Reliability of Estimator 

This time it is more difficult to calibrate the estimator Eq. (11) because our test models 
(Poisson and shuffled lattice) consist of Dirac 5's and zones of vanishing correlation. 

For the Poissonian case, £(r) = if r ^ 0, so that we expect £ es t(r) to fluctuate 
around zero with an amplitude proportional to the dispersion of the estimator, (£est( r ))- 
Absence of correlations allows an easy estimation of the dispersion: if N s is not too 
large (so that the probability that shells overlap is small), the numbers Tj in Eq. (12) 
are uncorrelated of each other and have a Poissonian distribution. One can then show 
immediately: 

(i\> = ( e ), (r^) - (r,)(r,) = (13) 



and then 
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Figure 5: Scaling relations for |£ e st(?") I when applied to a set of Poisson distributions and 
of shuffled lattices with spacing lh" 1 Mpc. The error bars for the Poisson distribution 
are la and for the lattice too small to be presented. The amplitudes are again arbitrary. 



We took dr oc r (logarithmic binning), so that we expect the amplitude of the fluctu- 
ations in £ es t(r) to decay like r~ 3 / 2 , as indeed observed in Fig. 5, where the error bars 
are again la errors when averaging over the ten random sets. The Figure also shows 
|£ est (r)| for the shuffled lattice with grid spacing of lh^ 1 Mpc. We believe that the 
observed r~ 2 -decay is again due to (£est( r ))> like in the Poissonian case. 



3.3.2 Application to iV-body Data 

Fig. 6 shows the result of applying the estimator Eq. (11) to the actual iV-body data at 
the initial redshift z = 50. We plot the absolute value as the correlation function 

tends to oscillate around zero, too. The curve is the average taken over the ten runs 
(as usual), but we do not plot the error bars as the data already show a noticeable 
level of noise. This noise is in fact so strong as to mask the signal (the ACDM behavior 
in this case); we already found this problem with the test models especially for the 
lattice distribution upon which the ICs are based on (cf. Eq. (1)). It seems that an 
improvement of the estimator (11) is required to extract reliable information in these 
extreme cases. 

Fig. 7 shows the same quantity for the z = data, where the analytical curves 
are the correlations derived from the linearly extrapolated P(k) (dotted line) and the 
non-linear P(k) from PD96 (solid line) in Eq. (10), respectively. This time we find 
a deviation of the estimated £(r) from the one predicted using PD96 on small scales. 
However, the error bars are 'only' la and the PD96 prediction still lies within the 3a 
level. One must also note that the estimator Eq. (11) is biased towards high density 
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Figure 6: Two-point correlation function for initial particle distribution at redshift 
z = 50. For clarity no error bars are shown due to a high level of noise. The solid line 
is the expected £(r) as given by Eq. (10). We plot |£(r)|, as the estimated function 
tends to oscillate around zero for scales smaller than the Nyquist wavelength (indicated 
by the vertical line). 

regions where most of the particles at z = will reside, since the shells are centered at 
(randomly chosen) particle positions rather than placing them randomly at any point 
(like the estimator Eq. (7) does, which also explains why the scatter for £est( r ) a ^ these 
small scales is much smaller than for o\j est (r) observed in Fig. 4). However, we varied 
the number of spheres N s from 50 to 100000 and could only detect a mild (if any) 
dependence of the amplitude on N s . 

Nevertheless, if one is to believe this discrepancy, then it is not obvious which one 
is to be blamed, the simulations or the PD96 fit. We have confirmed the excellent 
agreement of PD96 with our estimated P(k) in the probed range of wavenumbers 
(see Fig. 1). However, when using Eq. (10), one is extrapolating the PD96 fit to all 
wavenumbers. Clearly, the discrepancy should originate from the modes beyond the 
probed range, but with the information at hand, one cannot conclude whether their 
effect is estimated wrongly by the simulation or by the PD96 fit. 

3.4 Minkowski Functionals 

We have computed also the four scalar Minkowksi functionals (MF) of each realization 
(Mecke, Buchert & Wagner 1994). The MFs are morphological measures of the struc- 
ture, sensitive to correlations of order higher than the second. They include the genus 
statistics (Melott 1990) and have been used to quantify how filamentary or sponge-like 
the matter distribution looks like (Schmalzing et al. 1999), to study galaxy distribution 
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Figure 7: Same as Fig. 6 but for z — 0. The error bars are la errors when averaging 
over the ten runs. The two lines are again the analytical PD96 fit (solid) and the linear 
theory (dotted). 

in catalogs (Kerscher et al. 1997), and to address Gaussianity in the cosmic microwave 
background (Schmalzing & Gorski 1998). 

Like in deriving the power spectrum, we constructed a density field on a regular grid 
using the Triangular-Shape-Cloud method with two different resolutions: a 128 3 -cell 
grid and a 32 3 -cell grid. A density threshold was introduced and the boundary surface 
was determined between regions with a density below the threshold and regions with a 
density above it. Finally, the MFs of the boundary surface were determined. The four 
MFs are defined as follows: 

• M = volume enclosed by the surface, 

• Mi — area of the surface, 

• M 2 = integral over the surface of its mean curvature, 

• M 3 = integral over the surface of its Gaussian curvature, which coincides with 
the Euler characteristic (genus): 

M 3 = number of disconnected objects + number of holes — number of tunnels. 

Figure 8 shows the initial MFs as a function of the density threshold. The gen- 
eral shape of the plots can be explained qualitatively: for low thresholds, there is no 
bounding surface at all (due to the periodic boundary conditions), and hence M = 
total simulation volume, while M 1 = M 2 = M 3 = 0. As the threshold increases, there 
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Figure 8: MFs of the ten realizations of the initial matter distribution (z = 50) as a 
function of the density threshold 1 + 5, at two different spatial resolutions: 128 3 -cell 
grid (solid lines) and 32 3 -cell grid (dotted lines). 

appear first disconnected blobs of low density regions (M 3 > 0, increasing M 1 > 0, 
and the boundary surface is predominantly concave: M 2 < 0), later the blobs fuse 
together and tunnels arise (M 3 < 0), and finally the situation reverses and one ends 
up with independent blobs of high density regions (M 3 > 0, decreasing M 1 > 0, and 
the boundary surface is predominantly convex: M 2 > 0), until the threshold becomes 
larger than the maximum density (M = M 1 = M 2 = M 3 = 0). There are evidences 
that the zeros of M 3 (8) are strongly correlated with the percolation thresholds of the 
regions above or below the density threshold (Mecke & Wagner 1991). 

At the initial time, density fluctuations are small and Gaussian, which explains the 
symmetry of the MFs with respect to 5 = 0. However, a slight asymmetry can be 
detected for M 2 and M 3 on the 128 3 -cell grid: the spatial resolution is large enough 
that the MFs are sensitive to the finite-mass effects induced by the underlying point- 
particle distribution. Another difference between the two grids is the dispersion among 
realizations, which is larger when the spatial resolution is small. The scatter in the 
positions of the zeros and the values of the extrema tends to increase from M toward 
M 3 ; in the worst case, the value of the maxima of M 3 (5) has a 1-a error pa 10%. 

Figure 9 shows the MFs for the final matter distribution. The minimum density 
that can be resolved is 1 + 5 m i n ~ (mean interparticle distance/grid constant) 3 , which 
is 1 for the 128 3 -cell grid and 1/64 for the 32 3 -cell grid. Thus, the curves below these 
densities are in principle not reliable 5 . Apart from the asymmetry around 5 = 0, one 

5 For instance, the feature at 1 + 5 « 0.4 of M 2 and M 3 in the 128 3 -grid is likely a finite-mass effect 
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Figure 9: MFs of the ten realizations of the final matter distribution (z = 0) as a 
function of the logarithm of the density threshold log(l + 6), at two different spatial 
resolutions: 128 3 -cell grid (solid lines) and 32 3 -cell grid (dotted lines). 



observes in general that the scatter in the ordinate direction barely changes: the value 
of the maximum of M 3 (5) has an error 15%. The abscissa-dispersion, however, is 
larger than at the initial time: so, e.g., the zeros of M 3 (<5) in the 32 3 -grid have an error 
w 6%, while the uncertainty at the initial time is just ^ 0.3%. 



4 Analysis II: Dark Matter Halos 

The remaining analysis is going to focus on gravitationally bound halos, identified using 
the Bound-Density-Maxima method (BDM, Klypin & Holtzman 1997). We investigate 
the scatter in (large-scale) clustering patterns as well as internal properties of halos 
introduced by the random nature of the initial conditions. 



4.1 Identifying Halos 

We restricted our analysis to halos with more than 100 particles (M vir . imin > lO 12 /^ 1 M ). 
This lower mass limit can be used to derive a lower limit for the virial radius -R v ir,mm 
via 

47T 

M vir = yA ¥ir p fc ^ ir , (15) 

due to isolated particles (density peaks at a density w 1) in the voids of the structure. 
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Figure 10: Cumulative mass functions of BDM halos compared to Press-Schechter 
prediction (Press & Schechter 1974). The mass functions is the average taken over all 
ten runs and the error bars are la errors. 

where pb is the background density and A vir = 340 for the ACDM model under con- 
sideration. 

The BDM code identifies local overdensity peaks by smoothing the density field 
on a particular scale. The particle distribution was used to iteratively find potential 
halo centers as the centers of mass of 20,000 spheres with radius -R S pherc ~ -Rvir,min ~ 
162/i" 1 kpc centered about randomly chosen particles. Once the iteration converged 
for all spheres we repeated the procedure using successively smaller sphere radii down 
to 70/2T 1 kpc, about three times the force resolution. For each of these halo centers we 
stepped out in radial bins until the density drops below pbin < ^vhPb- This defined 
the outer radius _R v i r of the halo 6 . We discarded all halos with less than 100 particles 
within R vir for the further analysis. 7 

4.2 Mass Function of Halos 

The first quantity to investigate is the mass spectrum. We calculated the cumulative 
mass function n(>M) for our BDM halos and compared it to the analytical prediction 
of Press & Schechter (1974): 

6 If we want to identify halos-within-halos this method needs to be adjusted to account for the fact 
that the actual density of a satellite galaxy might not drop below A v ; r ph 

7 To cross check the completeness of our BDM halo catalogues we also performed a FOF analysis 
which shows a nearly 100% agreement and only an incompleteness in the BDM catalogues for halos 
less massive than 100 particles. 
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dM \ -k M a M 

where the variance ou is given by Eq. (6) and 5 C = 1.68. 

Fig. 10 shows that the average mass function of all ten runs is in good agree- 
ment with the PS prediction, what has been noted already by several other authors 
(Efstathiou et al. 1988; White, Efstathiou & Frenk 1993; Gross et al. 1998; Gover- 
nato et al. 1999; Jenkins et al. 2001). This again is another indicator that the initial 
conditions as well as the evolution by iV-body simulation are in fair agreement with 
theoretical predictions based on the analytical power spectrum and its time evolution. 
The discrepancy of the numerical n(> M) with the PS prediction at the low and high 
mass end of the mass function is also a well known fact (e.g. Governato et al. 1999) 
and not related to unreliable ICs or wrong iV-body modeling. Anyway, we are more 
interested in the scatter stemming from the random nature of the initial conditions. 
We are driven by the question, to what extent a single cosmological simulation can be 
representative for the volume under investigation. We observe that the scatter grad- 
ually increases from around 4% at the very low mass end resolved to about 50% for 
the most massive objects found in the simulation. According to the PS prediction, the 
scatter due to cosmic variance should enter via the amplitude predominantly, not 
the phases % of the ICs, Eq. (4). The observed increase of scatter with mass is then 
naturally explained also by the PS formula, given that larger masses are more sensitive 
to larger scales. 



dln<jM 
d\nM 



exp 



2< 



dM 
J AT 



(16) 



4.3 Halo-Halo Correlation Function 

The calculation of the halo-halo two-point correlation function is based on the estimator 
Eq. (11) again. However, this time we applied it only to the 500 most massive objects in 
the runs, which means fixing the number density of halos to nhaio = 2-10~ 3 (/i -1 Mpc)~ 3 . 
This choice for nhaio restricts the masses of the objects used from M ~ 3 • 10 14 /i _1 M 
down to M ~ 2 • lO 12 ^" 1 M . The result for the average taken over the ten BDM 
catalogs at redshift z = is shown in Fig. 11. The mean correlation function £ e st( r ) 
was fitted to a power law, 

£(r) = (r /ry , (17) 

over the range r G [0.5,20]/i -1 Mpc with the parameters r =4.26 ±0.44/i _1 Mpc and 
7=1.80 ±0.17. The la errors are of the order of 10% and indicate again only a mild 
dependence of the halo-halo correlation function on the variance introduced by the 
random nature of the initial conditions. Even though the scatter for the fundamental 
mode is ~ 20%, it does only marginally affect the statistical clustering properties of 
dark matter halos in the respective mass range. 
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Figure 11: Two-point correlation function £(r) for BDM halos at redshift z — 0. Error 
bars are again la. The thin solid line is the fit to a power law as given by Eq. (17). 




ioo iooo 

v [/r 1 kpc] 

Figure 12: Average density profile for the most massive BDM halo with la error bars 
along with a NFW and a Moore profile fit to the data. 
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4.4 Internal Properties of the Most Massive Halo 

Even though we are only resolving approximately 2% of the virial radius of the most 
massive particle group (cf. Table 1), we fitted a Navarro, Frenk & White (NFW) profile 
(Navarro, Frenk & White 1997), 

PNFw(r) K r/r.(l + r/r.)" ^ 
as well as a Moore profile (Moore et al. 1999), 

PMoorc(r) oc - 1 (19) 

(r/r s y- b (l + {r/r s ) Lb ) 

to the most massive halo found in the BDM catalogues. The question we are interested 
in is whether the scatter due to the random nature of the initial conditions can be 
made responsible for the difference in the central slope of the density profile described 
by those two fitting formulae. And from Fig. 12 we deduce that at least down to the 
resolved scale of 2% of the virial radius both analytical descriptions for the density 
profile do give indistinguishable good fits to the actual data; they both lie within the 
la error bars. However, we must stress that both profiles start to deviate even stronger 
from each other for even smaller scales not covered by the current study. Moreover, 
the reduced x 2 value for the NFW fit is marginally better than for the Moore fit, as 
one might anticipate from the behavior at small r values in Fig. 12. 

We conclude the analysis with Table 1 summarizing some internal properties calcu- 
lated for the most massive halo, i.e. mass M, circular velocity f C irc, velocity dispersion 
a v , virial radius r vir , concentration parameter 

c = r vir /r s , (20) 

where r s is the scale radius derived from the fit to the NFW profile Eq. (18), the spin 
parameter 

A = jJ\E\/(GM 5/2 ), (21) 

and the triaxiality parameter 



or — cr 



where a > b > c are the eigenvalues of the inertia tensor. 

This Table shows that the la variance for nearly all quantities is of the order of 
20%, like the variance of the fundamental mode k m i n = 2ir/B (cf. Fig. 1). Only the 
mass and the spin parameter show a larger scatter. 
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Table 1: Internal (averaged) properties of the most massive halo when averaged over 
ten runs. The errors are the la- value again. 

property variance 



M 


= (3.07 


± 1.60) 


lO 14 /^ 1 M 


52 % 


^circ 


= (1131 


± 199) 


km / sec 


18 % 


a v 


= (1172 


± 195) 


km/sec 


17% 




= (1344 


± 163) 


hr x kpc 


12 % 


C 


= 4.10 


± 0.91 




22 % 


A 


= 0.033 


± 0.018 




53 % 


T 


= 0.762 


± 0.102 




13 % 



5 Summary and Conclusions 

We presented the study of ten random realizations of a density field characterized by a 
cosmological power spectrum P(k) at redshift z = 50. These initial conditions for N- 
body simulations were tested with respect to their correlation properties. Recent claims 
by Baertschiger & Sylos Labini (2002) throw doubts on the ability of the commonly used 
method for generating the initial density field using particles (i.e. glass or grid preinitial 
distribution + the Zeldovich approximation, Eftstahiou et al. 1985) to clearly reproduce 
the analytical input correlations. The power spectrum P(k) and the mass variance 
(Tm(j') do not deviate from the expected behavior (including expected departures from 
the desired ACDM behavior due to finite mass and finite size effects). The estimated 
2-point correlation £(r) is too noisy to be used as a reliable credibility check; one 
cannot claim either that it reproduces the desired ACDM behavior or that it exhibits 
systematic deviations thereof. 

These initial conditions were then evolved forward in time until redshift z = 
using the publicly available adaptive mesh refinement code MLAPM (Knebe, Green & 
Binney 2001). This allowed us to explore the cosmic variance stemming from the ran- 
dom nature of the initial conditions, i.e., the scatter between different realizations of 
statistically identical initial conditions. We addressed the morphological properties of 
the matter distribution with the four Minkowski functionals as functions of a density 
threshold. The scatter grows in time, the one exhibiting a larger dispersion being the 
genus, of the order of 10% at z — 0. We also investigated the internal properties of 
DM halos, which have already been shown by other groups to be profoundly influenced 
by the surrounding large-scale structure, which in turn is sensitive to A;-modes ~ fun- 
damental mode (Colberg et al. 1999). We find that the scatter in the properties of the 
most massive object(s) forming in the box is ~ 20%, and as high as ~ 50% for some 
properties such as the mass or the spin parameter. 

An interesting question is whether this scatter is induced mainly by the cosmic 
variance of the amplitude at scales around the fundamental mode, or by the cosmic 
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variance of the random phases. There is certainly a propagation of the error in the 
initial large-scale amplitude by power transfer towards smaller scales. In fact, when 
comparing our data to the (non-)linear fit of Peacock & Dodds (1996) for the power 
spectrum and to the prediction by Press & Schechter (1974) for the mass function, 
we find good agreement. The derivation of both results is based on the hypothesis of 
small influence from coupling of modes at some k to modes with larger k; our results 
support this assumption, as far as the statistical estimators we probed are concerned. It 
would now be interesting to investigate in detail the actual influence of the large waves 
onto the small scale structure. This would also shed some light on the credibility of 
running small simulation boxes to very low redshifts as already done by several groups 
(e.g. Dave et al. 2001, Avila-Reese et al. 2001, Gnedin 2000, Colin, Avila-Reese & 
Valenzuela 2000), but we leave this to a future study. 
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